#include <stdlib.h>
#include <petsc.h>
#include "vtkDoubleArray.h"
#include "vtkPointData.h"
#include "vtkUnstructuredGrid.h"
#include "vtkXMLUnstructuredGridWriter.h"
#include "par_distribution.h"
#include "stateproblem.h"
#include "matrix_umf.h"
#include "output.h"
#include "integration/quadrature.h"
#include "integration/constpoint.h"
#include "integration/linpoint.h"
#include "integration/linbubblepoint.h"
#include "integration/lindiscpoint.h"
#include "integration/nclinpoint.h"
#include "integration/quadpoint.h"
#include "integration/quaddiscpoint.h"
#include "integration/quadbubblepoint.h"
#include "fem.h"


using namespace mu;


void start_clock(clock_t &t)
{
  t = clock();
}

double stop_clock(clock_t &t)
{
  return ((double)(clock()-t))/CLOCKS_PER_SEC;
}

string replace(string s, map<string,string> &subs)
{
  size_t found;
  
  for (map<string,string>::iterator it=subs.begin(); it!=subs.end(); it++)
  {
    found = s.find(it->first);
    while (found != string::npos)
    {
      s.replace(found, it->first.length(), it->second);
      found = s.find(it->first);
    }
  }
  
  return s;
}

CComponent::CComponent()
{
}

CComponent::~CComponent()
{
  for (map<string,Parser*>::iterator pit=p.begin(); pit!=p.end(); pit++)
  {
    delete pit->second;
  }
  p.clear();
    
  delete[] node_value;
  delete[] node_grad;
  delete[] node_count;
}

CUnknown::CUnknown()
{
}

CUnknown::~CUnknown()
{
  for (map<string,Parser*>::iterator pit=normalbc.begin(); pit!=normalbc.end(); pit++) delete pit->second;
  normalbc.clear();
  for (map<string,Parser*>::iterator pit=tangentbc.begin(); pit!=tangentbc.end(); pit++) delete pit->second;
  tangentbc.clear();
}


PPdata::PPdata(string _name, int _type, vector<string> _expr, int size) :
  name(_name),
  type(_type),
  expressions(_expr)
{
  for (int i=0; i<type; i++)
    data.push_back(new double[size]);
}

PPdata::~PPdata()
{
  for (int i=0; i<type; i++)
    delete[] data[i];
}

string &PPdata::getExpr(int i)
{
  return expressions[i];
}

double *PPdata::getData(int i)
{
  return data[i];
}


CStateProblem::CStateProblem() :
  mesh(0),
  Q3d(0),
  Q2d(0),
  wm(0),
  rstop(1e-6),
  tau(0),
  tmax(0),
  dim(3),
  loaded_solution(0),
  calculated_solution(0),
  logfile(0),
  output(0)
{
  IP3D.clear();
  IP2D.clear();
  IP1D.clear();
  U.clear();
  Uinfo.clear();
  PPfields.clear();
  
  wm = new CPlainWManager;
}


CStateProblem::~CStateProblem()
{
  for (Type<3>::t_IP::iterator it=IP3D.begin(); it != IP3D.end(); it++) delete (*it).second;
  IP3D.clear();
  for (Type<2>::t_IP::iterator it=IP2D.begin(); it != IP2D.end(); it++) delete (*it).second;
  IP2D.clear();
  for (Type<1>::t_IP::iterator it=IP1D.begin(); it != IP1D.end(); it++) delete (*it).second;
  IP1D.clear();
  
  U.clear();
  
  for (int i=0; i<(int)Uinfo.size(); i++) delete Uinfo[i];
  Uinfo.clear();
  
  if (mesh == 0)
  {
    delete mesh;
    mesh = 0;
  }
  
  for (int i=0; i<q.size(); i++) drq[i].clear();
  delete[] drq;
  
  delete[] r;
  
  delete[] f_array;

  delete[] petsc_dof_numbering;

  params.clear();
  
  for (map<string,PPdata *>::iterator it = PPfields.begin(); it != PPfields.end(); it++) delete it->second;
  PPfields.clear();
  
  if (wm != 0) delete wm;
  
  if (output != 0) delete output;
  
  if (logfile != 0) fclose(logfile);
}

void CStateProblem::getSubstitutions(const char *argument)
{
  unsigned int pos = 0, pos_equality;
  string s = argument, key, data;
  
  while (pos < s.length())
  {
    pos_equality = s.find_first_of('=', pos);
    key = "%" + s.substr(pos, pos_equality-pos);
    pos = s.find_first_of(' ', pos_equality);
    data = s.substr(pos_equality+1,pos-pos_equality-1);
    substitutions[key] = data;
    pos = s.find_first_not_of(' ', pos);
  }
}

void CStateProblem::readConfig(char *file)
{
  ifstream f(file);
  string l, key, data, data1, data2;
  
  // check if file exists
  if ((bool)f == false) return;
  
  log(W_PROG, "Reading config file %s.\n", file);
  
  while (!f.eof())
  {
    getline(f, l);
    key = l.substr(0, l.find_first_of('='));
    data = replace(l.substr(l.find_first_of('=')+1), substitutions);
    
    if (key.compare("") != 0 && key.compare(0, 1, "#") != 0)
    { // unrecognized keys are stored
      params[key] = data;
    }
  }
  
  for (map<string,string>::iterator it=params.begin(); it!=params.end(); it++)
  {
    key = it->first;
    data = it->second;
    
    if (key.compare("meshfile") == 0)
    {
      int start_pos = data.find_first_of(',')+1;
      unsigned int comma_pos = data.find_first_of(',', start_pos);
      data1 = data.substr(0, start_pos-1);
      data2 = data.substr(start_pos, comma_pos-start_pos);
      string s_cells;
      SaveCells cells = SaveNone;
      do
      {
    	  start_pos = comma_pos+1;
    	  comma_pos = data.find_first_of(',', start_pos);
    	  s_cells = data.substr(start_pos, comma_pos-start_pos);
    	  if (s_cells.compare("3d") == 0)
    	  {
    		  cells |= Save3D;
    	  }
    	  else if (s_cells.compare("2d") == 0)
    	  {
    		  cells |= Save2D;
    	  }
    	  else if (s_cells.compare("1d") == 0)
    	  {
    		  cells |= Save1D;
    	  }
      } while (comma_pos < data.size());

      if (data2.compare("linear") == 0)
      {
        loadMesh(FE_LIN, data1.c_str(), cells);
      }
      else if (data2.compare("quadratic") == 0)
      {
        loadMesh(FE_QUAD, data1.c_str(), cells);
      }
      else
      {
        log(W_PROG, "Error reading config file %s: Unsupported mesh type (%s).\n", file, data2.c_str());
	throw 1;
	exit(1);
      }
    }
    else if (key.compare("logfile") == 0)
    { // open log file
      logfile = fopen(data.c_str(), "w");
    }
    else if (key.compare("rstop") == 0)
    { // stopping tolerance for residual norm
      rstop = atof(data.c_str());
    }
    else if (key.compare("timestep") == 0)
    { // open log file
      tau = atof(data.c_str());
    }
    else if (key.compare("stoptime") == 0)
    { // open log file
      tmax = atof(data.c_str());
    }
    else if (key.compare(0, 11, "save_scalar") == 0)
    { // save post process field
      string name = key.substr(12, key.size()-13);
      addPPscalar(name, data);
    }
  }
  afterReadConfig();
}

void CStateProblem::loadMesh(int meshtype, const char *file, SaveCells cells)
{
	  int ierr, my_proc, num_of_procs;
	  ierr=MPI_Comm_rank(PETSC_COMM_WORLD, &(my_proc));
	  ASSERT( ! ierr , "Can not get MPI rank.\n" );
	  ierr=MPI_Comm_size(PETSC_COMM_WORLD, &(num_of_procs));
	  ASSERT( ! ierr  , "Can not get MPI size.\n" );

	  switch (meshtype)
	  {
	  case FE_LIN:
		  mesh = new LinMesh;
		  break;
	  case FE_QUAD:
		  mesh = new QuadMesh;
		  break;
	  default:
		  log(W_PROG, "Error: Unknown mesh type: %d.\n", meshtype);
		  throw 1;
		  exit(1);
	  }
	  mesh->readGMSH(file, (num_of_procs==1)?0:my_proc+1);
  
	  afterLoadMesh();
  
	  log(W_PROG, "Loaded mesh file %s:\n %d volumes,\n %d faces,\n %d edges,\n %d points,\n %d nodes,\n %d nodes shared by partitions.\n",
			  file,
			  mesh->getNvolumes(),
			  mesh->getNfaces(),
			  mesh->getNedges(),
			  mesh->getNpoints(),
			  mesh->getNnodes(),
			  mesh->getSharedNodes().size());

	  // Initialize output writer class
	  output = new Output(mesh, cells);
}

void CStateProblem::setDimension(int dimension)
{
  dim = dimension;
}



void CStateProblem::registerQuadratures(const Quadrature<3> *_q3d, const Quadrature<2> *_q2d, const Quadrature<1> *_q1d)
{
  Q3d = _q3d;
  Q2d = _q2d;
  Q1d = _q1d;
}

void CStateProblem::addComponent(int utype, bool _set_mv, int _mv_segment, double _mv)
{
  CComponent *u = new CComponent;
  Parser *p;
  char no[5];
  
  switch (utype)
  {
    case FE_CONST:
      if (IP3D.find(FE_CONST) == IP3D.end())
      {
	IP3D[FE_CONST] = new ConstPoint<3>(mesh, Q3d);
	IP2D[FE_CONST] = new ConstPoint<2>(mesh, Q2d);
	IP1D[FE_CONST] = new ConstPoint<1>(mesh, Q1d);
      }
      break;
    case FE_LIN:
      if (IP3D.find(FE_LIN) == IP3D.end())
      {
        if (dynamic_cast<LinMesh *>(mesh) != 0)
        {
	  IP3D[FE_LIN] = new LinPoint<3>(mesh, Q3d);
	  IP2D[FE_LIN] = new LinPoint<2>(mesh, Q2d);
	  IP1D[FE_LIN] = new LinPoint<1>(mesh, Q1d);
	}
	else if (dynamic_cast<QuadMesh *>(mesh) != 0)
	{
	  IP3D[FE_LIN] = new LinPointOnQuadMesh<3>(mesh, Q3d);
	  IP2D[FE_LIN] = new LinPointOnQuadMesh<2>(mesh, Q2d);
	  IP1D[FE_LIN] = new LinPointOnQuadMesh<1>(mesh, Q1d);
	}
	else
	{
	  log(W_PROG, "Error: Component type is incompatible with mesh.\n");
	  throw 1;
	  exit(1);
	}
      }
      break;
    case FE_LIN_DISC:
      if (IP3D.find(FE_LIN_DISC) == IP3D.end())
      {
        IP3D[FE_LIN_DISC] = new LinDiscPoint<3>(mesh, Q3d);
        IP2D[FE_LIN_DISC] = new LinDiscPoint<2>(mesh, Q2d);
        IP1D[FE_LIN_DISC] = new LinDiscPoint<1>(mesh, Q1d);
      }
      break;
    case FE_LIN_NC:
      if (IP3D.find(FE_LIN_NC) == IP3D.end())
      {
        IP3D[FE_LIN_NC] = new NCLinPoint<3>(mesh, Q3d);
        IP2D[FE_LIN_NC] = new NCLinPoint<2>(mesh, Q2d);
        IP1D[FE_LIN_NC] = new LinPoint<1>(mesh, Q1d);
      }
      break;
    case FE_LIN_BUBBLE:
      if (IP3D.find(FE_LIN_BUBBLE) == IP3D.end())
      {
        IP3D[FE_LIN_BUBBLE] = new LinBubblePoint<3>(mesh, Q3d);
        IP2D[FE_LIN_BUBBLE] = new LinBubblePoint<2>(mesh, Q2d);
		IP1D[FE_LIN_BUBBLE] = new LinBubblePoint<1>(mesh, Q1d);
      }
      break;
    case FE_QUAD:
      if (IP3D.find(FE_QUAD) == IP3D.end())
      {
        IP3D[FE_QUAD] = new QuadPoint<3>(mesh, Q3d);
        IP2D[FE_QUAD] = new QuadPoint<2>(mesh, Q2d);
		IP1D[FE_QUAD] = new QuadPoint<1>(mesh, Q1d);
      }
      break;
    case FE_QUAD_DISC:
      if (IP3D.find(FE_QUAD_DISC) == IP3D.end())
      {
        IP3D[FE_QUAD_DISC] = new QuadDiscPoint<3>(mesh, Q3d);
        IP2D[FE_QUAD_DISC] = new QuadDiscPoint<2>(mesh, Q2d);
        IP1D[FE_QUAD_DISC] = new QuadDiscPoint<1>(mesh, Q1d);
      }
      break;
    case FE_QUAD_BUBBLE:
      if (IP3D.find(FE_QUAD_BUBBLE) == IP3D.end())
      {
        IP3D[FE_QUAD_BUBBLE] = new QuadBubblePoint<3>(mesh, Q3d);
        IP2D[FE_QUAD_BUBBLE] = new QuadBubblePoint<2>(mesh, Q2d);
        IP1D[FE_QUAD_BUBBLE] = new QuadPoint<1>(mesh, Q1d);
      }
      break;
    default:
      log(W_PROG, "Error: Component type %d not recognized.\n", utype);
      throw 1;
      exit(1);
  }
  u->ip3d = IP3D[utype];
  u->ip2d = IP2D[utype];
  u->ip1d = IP1D[utype];
  
  u->set_mean_value = _set_mv;
  if (_set_mv)
  {
    u->mv_segment = _mv_segment;
    u->mean_value = _mv;
    u->volume     = 0;
    u->mean = new double[u->ip3d->getMeshNdof()];
    for (int j=0; j<u->ip3d->getMeshNdof(); j++) u->mean[j] = 0;
  }
  
  string bctag;
  for (map<int,string>::iterator bcname=mesh->beginPhysicalNames(); bcname!=mesh->endPhysicalNames(); bcname++)
  {
    sprintf(no, "%d", (int)U.size());
    bctag = "bc[" + bcname->second + "," + no + "]";
    if (params.find(bctag) != params.end())
    {
      p = new Parser;
      p->DefineVar("x", &x_d);
      p->DefineVar("y", &y_d);
      p->DefineVar("z", &z_d);
      p->DefineVar("t", &time);
      p->SetExpr(params[bctag]);
      u->p[bcname->second] = p;
      params.erase(bctag);
    }
  }
  
  int N = output->size();
  
  u->node_value = new double[N];
  u->node_grad  = new double[3*N];
  u->node_count = 0;
  
  U.push_back(u);
}

void CStateProblem::addScalar(string name, int utype, bool _set_mv, int _mv_segment, double _mv)
{
  CUnknown *ui = new CUnknown;
  ui->name  = name;
  ui->type  = UT_SCALAR;
  ui->index = U.size();
  
  Uinfo.push_back(ui);
  
  addComponent(utype, _set_mv, _mv_segment, _mv);
  int d = output->registerNodeData(name, UT_SCALAR);
  output->addComponent(d, U[ui->index]->node_value);
}

void CStateProblem::addVector(string name, int utype, bool _set_mv, int _mv_segment, double _mv)
{
  CUnknown *ui = new CUnknown;
  ui->name  = name;
  ui->type  = UT_VECTOR;
  ui->index = U.size();
  
  Uinfo.push_back(ui);
  
  int d = output->registerNodeData(name, UT_VECTOR);
  
  for (int i=0; i<dim; i++)
  {
    addComponent(utype, _set_mv, _mv_segment, _mv);
    output->addComponent(d, U[ui->index+i]->node_value);
  }
  
  string bctag;
  Parser *p;
  for (map<int,string>::iterator bcname=mesh->beginPhysicalNames(); bcname!=mesh->endPhysicalNames(); bcname++)
  {
    bctag = "normalbc[" + name + "," + bcname->second + "]";
    if (params.find(bctag) != params.end())
    {
      p = new Parser;
      p->DefineVar("x", &x_d);
      p->DefineVar("y", &y_d);
      p->DefineVar("z", &z_d);
      p->DefineVar("t", &time);
      p->SetExpr(params[bctag]);
      ui->normalbc[bcname->second] = p;
      params.erase(bctag);
    }
    bctag = "tangentbc[" + name + "," + bcname->second + "]";
    if (params.find(bctag) != params.end())
    {
      p = new Parser;
      p->DefineVar("x", &x_d);
      p->DefineVar("y", &y_d);
      p->DefineVar("z", &z_d);
      p->DefineVar("t", &time);
      p->SetExpr(params[bctag]);
      ui->tangentbc[bcname->second] = p;
      params.erase(bctag);
    }
  }
}

void CStateProblem::addSymTensor(string name, int utype, bool _set_mv, int _mv_segment, double _mv)
{
  CUnknown *ui = new CUnknown;
  ui->name  = name;
  ui->type  = UT_SYM_TENSOR;
  ui->index = U.size();
  
  Uinfo.push_back(ui);
  
  int d = output->registerNodeData(name, UT_SYM_TENSOR);
  
  int n = 0;
  switch (dim)
  {
    case 2:
      n = 3;
      break;
    case 3:
      n = 6;
      break;
  }
  for (int i=0; i<n; i++)
  {
    addComponent(utype, _set_mv, _mv_segment, _mv);
    output->addComponent(d, U[ui->index+i]->node_value);
  }
}

void CStateProblem::addTensor(string name, int utype, bool _set_mv, int _mv_segment, double _mv)
{
  CUnknown *ui = new CUnknown;
  ui->name  = name;
  ui->type  = UT_TENSOR;
  ui->index = U.size();
  
  Uinfo.push_back(ui);
  
  int d = output->registerNodeData(name, UT_TENSOR);
  
  for (int i=0; i<dim*dim; i++)
  {
    addComponent(utype, _set_mv, _mv_segment, _mv);
    output->addComponent(d, U[ui->index+i]->node_value);
  }
}

int CStateProblem::getProblemNdof()
{
  int n = 0;
  for (int i=0; i<(int)U.size(); i++)
  {
    n += U[i]->ip3d->getMeshNdof();
  }
  return n;
}

void CStateProblem::prepare()
{
    int ierr, my_proc, num_of_procs;
    ierr=MPI_Comm_rank(PETSC_COMM_WORLD, &(my_proc));
    ASSERT( ! ierr , "Can not get MPI rank.\n" );
    ierr=MPI_Comm_size(PETSC_COMM_WORLD, &(num_of_procs));
    ASSERT( ! ierr  , "Can not get MPI size.\n" );

    q.resize(getProblemNdof());
    drq = new t_drq[q.size()];
    r   = new double[q.size()];
  
    // Create parallel distribution of dofs
    if (num_of_procs > 1)
    {
    	double lsize = 0;
    	for (unsigned int i=0; i<U.size(); i++)
    	{
    		for (int j=0; j<U[i]->ip3d->getMeshNdof(); j++)
    		{
    			if (mesh->getNodePartNo(j) == my_proc+1)
    			{
    				lsize++;
    			}
    		}
    	}
    	dof_distribution = new Distribution(lsize);
    }
    else
    {
    	dof_distribution = new Distribution(q.size());
    }

    // create dof numbering in the PETSc matrix
    petsc_dof_numbering = new int[q.size()];
    int j=0;
    int poffset = (dof_distribution->np() == 1)?0:1;
    int offset;
    for (int p=0; p<dof_distribution->np(); p++)
    {
    	offset = 0;
    	for (unsigned int c=0; c<U.size(); c++)
    	{
    		for (set<int>::iterator i=mesh->getPartitionNodes(p+poffset).begin(); i!=mesh->getPartitionNodes(p+poffset).end(); i++)
    		{
    			petsc_dof_numbering[offset+*i] = j;
    			j++;
    		}
    		offset += U[c]->ip3d->getMeshNdof();
    	}
    }

    // create local vector of rhs/solution
    f_array = new double[dof_distribution->lsize()];
  
    log(W_PROG, "Allocated data structures (problem size %d, local sizes ", q.size());
    for (int i=0; i<dof_distribution->np()-1; i++)
    	log(W_PROG, "%d, ", dof_distribution->lsize(i));
    log(W_PROG, "%d).\n", dof_distribution->lsize(dof_distribution->np()-1));

    // load solution
    if (loadSolution())
    {
    	log(W_PROG, "Loaded solution from file %s.\n", params["loadfile"].c_str());
    }
    else
    {
    	q = 0;
    	log(W_PROG, "Starting from initial guess.\n");
    }
    q.setIndependent();
  
    afterPrepare();
}

void CStateProblem::addResidual(CVar &r_loc, int ind)
{
  r[ind] += r_loc.getVal();
  for (int j=0; j<r_loc.getNder(); j++)
  {
    if (r_loc.getDerVar(j) >= q(0).getInd() && r_loc.getDerVar(j) < q(0).getInd()+q.size())
    {
      drq[ind][r_loc.getDerVar(j)-q(0).getInd()] += r_loc.getDer(j);
    }
  }
}

void CStateProblem::setResidual(CVar &r_loc, int ind)
{
  r[ind] = r_loc.getVal();
  drq[ind].clear();
  for (int j=0; j<r_loc.getNder(); j++)
  {
    if (r_loc.getDerVar(j) >= q(0).getInd() && r_loc.getDerVar(j) < q(0).getInd()+q.size())
    {
      drq[ind][r_loc.getDerVar(j)-q(0).getInd()] = r_loc.getDer(j);
    }
  }
}

void CStateProblem::clearResidual(int ind)
{
  r[ind] = 0;
  drq[ind].clear();
}

void calculateTangentsFromNormal(CVec &n, CVec &t1, CVec &t2)
{
  CVar norm;
  t1.resize(3);
  t2.resize(3);
  if (fabs(n(2).getVal()) > 1e-3)
  {
    norm = sqrt(n(0)*n(0) + n(2)*n(2))/n(2);
    t1(0) = (CVar)1/norm;
    t1(1) = 0;
    t1(2) = -n(0)/(n(2)*norm);
  }
  else if (fabs(n(1).getVal()) > 1e-3)
  {
    norm = sqrt(n(0)*n(0) + n(1)*n(1))/n(1);
    t1(0) = (CVar)1/norm;
    t1(1) = -n(0)/(n(1)*norm);
    t1(2) = 0;
  }
  else
  {
    norm = sqrt(n(0)*n(0) + n(2)*n(2))/n(0);
    t1(0) = -n(2)/(n(0)*norm);
    t1(1) = 0;
    t1(2) = (CVar)1/norm;
  }
  t2(0) = n(1)*t1(2) - n(2)*t1(1);
  t2(1) = n(2)*t1(0) - n(0)*t1(2);
  t2(2) = n(0)*t1(1) - n(1)*t1(0);
}

void CStateProblem::assembleVolumeIntegral()
{
	int c, e, i, ind, k, offset, pct=0;
	IntPoint<3> *ip;
	CVar r_loc;
	Type<3>::t_IP::iterator it;

	log(W_SOLVE, "-Assembling volume integrals: ");

    // assemble volume integrals
    for (e=0; e<mesh->getNvolumes(); e++)
    {
    	// skip volumes of different partitions
    	if (dof_distribution->np() > 1 && mesh->getVolumePartNo(e) != dof_distribution->myp()+1) continue;

      // calculate integration point data
      for (it=IP3D.begin(); it!=IP3D.end(); it++) it->second->calcE(e, EVAL_FG);
      for (k=0; k<Q3d->n; k++)
      {
        // calculate integration point data
        for (it=IP3D.begin(); it!=IP3D.end(); it++) it->second->calcQP(k, EVAL_FG);
        // evaluate unknowns at integration points
        offset = 0;
        for (c=0; c<(int)U.size(); c++)
        {
		  U[c]->ip3d->eval(q, e, offset, U[c]->val, U[c]->grad, EVAL_FG);
		  if (tau>0)
		  {
			U[c]->ip3d->eval(oq, e, offset, U[c]->o_val, U[c]->o_grad, EVAL_FG);
			U[c]->ip3d->eval(ooq, e, offset, U[c]->oo_val, U[c]->oo_grad, EVAL_FG);
		  }
		  offset += U[c]->ip3d->getMeshNdof();
        }
		x = U[0]->ip3d->getX();
		y = U[0]->ip3d->getY();
		z = U[0]->ip3d->getZ();
		afterEvaluateVolumeComponents(e, k);

        // assemble element integrals
        for (it=IP3D.begin(); it!=IP3D.end(); it++)
        {
          ip = it->second;
          for (i=0; i<ip->getNdof(); i++)
          {
            ip->evalBasis(i, phi, gphi, EVAL_FG);
            ind = ip->mapDof(e, i);
			offset = 0;
			for (c=0; c<(int)U.size(); c++)
			{
			  if (U[c]->ip3d == ip && q(offset+ind).getDerOfVar(q(offset+ind).getInd()) == 1)
			  {
				r_loc = volumeIntegral(c, e, k)*ip->dx;
				addResidual(r_loc, ind+offset);
				if (U[c]->set_mean_value && U[c]->mv_segment == mesh->getVolumePhysicalNo(e))
				{
				  U[c]->mean[ind] += phi.getVal()*ip->dx;
				  U[c]->volume    += ip->dx/ip->getNdof();
				}
			  }
			  offset += U[c]->ip3d->getMeshNdof();
			}
		  }
        }
      }
      if (pct < (e*100/mesh->getNvolumes()))
      {
        pct = (e*100/mesh->getNvolumes());
        logPct(W_SOLVE, pct);
      }
    } // for e
    logPct(W_SOLVE, 100);
}

void CStateProblem::assembleFaceIntegral()
{
  log(W_SOLVE, "Assembling face integrals: ");
  int c, e, i, ind, indd, j, k, offset, pct=0;
  IntPoint<2> *ip;
  CVar r_loc;
  Type<2>::t_IP::iterator it;
    // assemble face integrals
    for (e=0; e<mesh->getNfaces(); e++)
    {
    	// skip volumes of different partitions
    	if (dof_distribution->np() > 1 && mesh->getFacePartNo(e) != dof_distribution->myp()+1) continue;

      // calculate integration point data
      for (it=IP2D.begin(); it!=IP2D.end(); it++) it->second->calcE(e, EVAL_FG);
      for (k=0; k<Q2d->n; k++)
      {
        offset = 0;
        // calculate integration point data
        for (it=IP2D.begin(); it!=IP2D.end(); it++) it->second->calcQP(k, EVAL_FG);
        // evaluate unknowns at integration points
        for (c=0; c<(int)U.size(); c++)
        {
          U[c]->ip2d->eval(q, e, offset, U[c]->val, U[c]->grad, EVAL_FG);
	  if (tau>0)
	  {
	    U[c]->ip2d->eval(oq, e, offset, U[c]->o_val, U[c]->o_grad, EVAL_FG);
	    U[c]->ip2d->eval(ooq, e, offset, U[c]->oo_val, U[c]->oo_grad, EVAL_FG);
	  }
	  if (U[c]->ip2d->isDiscontinuous())
	  {
	    for (i=0; i<mesh->getFaceNparents(e); i++)
              U[c]->ip2d->evalDG(q, e, i, offset, U[c]->dg_val[i], U[c]->dg_grad[i], EVAL_FG);
	  }
          offset += U[c]->ip2d->getMeshNdof();
        }
	x = U[0]->ip2d->getX();
  	y = U[0]->ip2d->getY();
	z = U[0]->ip2d->getZ();
	afterEvaluateFaceComponents(e, k);
	
        // assemble element integrals
        for (it=IP2D.begin(); it!=IP2D.end(); it++)
        {
          ip   = it->second;
          
          for (i=0; i<ip->getNdof(); i++)
          {
            ip->evalBasis(i, phi, gphi, EVAL_FG);
            ind = ip->mapDof(e, i);
	    offset = 0;
  	    for (c=0; c<(int)U.size(); c++)
    	    {
  	      if (U[c]->ip2d == ip && q(offset+ind).getDerOfVar(q(offset+ind).getInd()) == 1)
	      {
	        r_loc = faceIntegral(c, e)*ip->dx;
		addResidual(r_loc, ind+offset);
		if (U[c]->set_mean_value && U[c]->mv_segment == mesh->getFacePhysicalNo(e))
		{
		  U[c]->mean[ind] += phi.getVal()*ip->dx;
		  U[c]->volume    += ip->dx/ip->getNdof();
		}
	      }
	      if (U[c]->ip2d->isDiscontinuous() && U[c]->ip2d == ip)
	      {
	        for (j=0; j<mesh->getFaceNparents(e); j++)
	        {
	          indd = ip->mapDof(e, i, j);
	          if (q(offset+indd).getDerOfVar(q(offset+indd).getInd()) == 1)
	          {
	            r_loc = faceIntegral(c, e, j)*ip->dx;
	            addResidual(r_loc, indd+offset);
	          }
	        }
	      }
	      offset += U[c]->ip2d->getMeshNdof();
  	    }
	  }
        }
      }
      if (pct < (e*100/mesh->getNfaces()))
      {
        pct = (e*100/mesh->getNfaces());
        logPct(W_SOLVE, pct);
      }
    }
    logPct(W_SOLVE, 100);
}

void CStateProblem::assembleEdgeIntegral()
{
  int c, e, i, ind, indd, j, k, offset, pct=0;
  IntPoint<1> *ip;
  CVar r_loc;
  Type<1>::t_IP::iterator it;
  
  // assemble edge integrals
  log(W_SOLVE, "Assembling edge integrals: ");
  for (e=0; e<mesh->getNedges(); e++)
  {
  	// skip volumes of different partitions
  	if (dof_distribution->np() > 1 && mesh->getEdgePartNo(e) != dof_distribution->myp()+1) continue;

    // calculate integration point data
    for (it=IP1D.begin(); it!=IP1D.end(); it++) it->second->calcE(e, EVAL_FG);
    for (k=0; k<Q1d->n; k++)
    {
      offset = 0;
      // calculate integration point data
      for (it=IP1D.begin(); it!=IP1D.end(); it++) it->second->calcQP(k, EVAL_FG);
      // evaluate unknowns at integration points
      for (c=0; c<(int)U.size(); c++)
      {
    	U[c]->ip1d->eval(q, e, offset, U[c]->val, U[c]->grad, EVAL_FG);
	if (tau>0)
	{
	  U[c]->ip1d->eval(oq, e, offset, U[c]->o_val, U[c]->o_grad, EVAL_FG);
	  U[c]->ip1d->eval(ooq, e, offset, U[c]->oo_val, U[c]->oo_grad, EVAL_FG);
	}
	if (U[c]->ip1d->isDiscontinuous())
	{
	  for (i=0; i<mesh->getEdgeNparents(e); i++)
            U[c]->ip1d->evalDG(q, e, i, offset, U[c]->dg_val[i], U[c]->dg_grad[i], EVAL_FG);
	}
	offset += U[c]->ip1d->getMeshNdof();
      }
      x = U[0]->ip1d->getX();
      y = U[0]->ip1d->getY();
      z = U[0]->ip1d->getZ();
      afterEvaluateEdgeComponents(e, k);
	
      // assemble element integrals
      for (it=IP1D.begin(); it!=IP1D.end(); it++)
      {
        ip   = it->second;
        for (i=0; i<ip->getNdof(); i++)
        {
          ip->evalBasis(i, phi, gphi, EVAL_FG);
          
          ind = ip->mapDof(e, i);
	  offset = 0;
  	  for (c=0; c<(int)U.size(); c++)
    	  {
  	    if (U[c]->ip1d == ip && q(offset+ind).getDerOfVar(q(offset+ind).getInd()) == 1)
	    {
	      r_loc = edgeIntegral(c, e)*ip->dx;
	      addResidual(r_loc, ind+offset);
	      if (U[c]->set_mean_value && U[c]->mv_segment == mesh->getEdgePhysicalNo(e))
	      {
		U[c]->mean[ind] += phi.getVal()*ip->dx;
		U[c]->volume    += ip->dx/ip->getNdof();
	      }
	    }
	    if (U[c]->ip1d->isDiscontinuous() && U[c]->ip1d == ip)
	    {
	      for (j=0; j<mesh->getEdgeNparents(e); j++)
	      {
	        indd = ip->mapDof(e, i, j);
	        if (q(offset+indd).getDerOfVar(q(offset+indd).getInd()) == 1)
	        {
	          r_loc = edgeIntegral(c, e, j)*ip->dx;
	          addResidual(r_loc, indd+offset);
	        }
	      }
            }
	    offset += U[c]->ip1d->getMeshNdof();
	  }
	  
	}
      }
    }
    if (pct < (e*100/mesh->getNedges()))
    {
      pct = (e*100/mesh->getNedges());
      logPct(W_SOLVE, pct);
    }
  }
  logPct(W_SOLVE, 100);
}

void CStateProblem::assembleFaceNormalTangentBC()
{
  int offset, c, e, i, j, ind;
  IntPoint<2> *ip;
  
  if (dim != 3) return;
  
  for (map<int,string>::iterator b=mesh->beginPhysicalNames(); b != mesh->endPhysicalNames(); b++)
  {
    offset = 0;
    for (c=0; c<(int)Uinfo.size(); c++)
    {
      ip = U[Uinfo[c]->index]->ip2d;
      int udim=ip->getMeshNdof();
      if (Uinfo[c]->normalbc.find(b->second) != Uinfo[c]->normalbc.end() || Uinfo[c]->tangentbc.find(b->second) != Uinfo[c]->tangentbc.end())
      {
        map<int,int> index;
        map<int,CVec> n;
        CVec t1, t2, ln, lt1, lt2;
        index.clear();
        n.clear();
        for (e=0; e<mesh->getNfaces(); e++)
	{
	  if (mesh->getPhysicalName(mesh->getFacePhysicalNo(e)) == b->second)
	  {
	    mesh->calculateFaceNormal(e, ln, lt1, lt2);
	    for (i=0; i<ip->getNdof(); i++)
	    {
	      ind = ip->mapDof(e,i);
	      if (index.find(ind) == index.end())
	      {
	        index[ind] = mesh->getFaceNode(e, i);
		n    [ind] = ln;
	      }
	      else
	      {
		for (j=0; j<3; j++) n[ind](j) += ln(j);
	      }
	    }
	  }
	}
	CVar r1, r2, r3;
	double norm;
	t_drq d1, d2;
	int i1, i2, i3;
        for (map<int,int>::iterator it=index.begin(); it != index.end(); it++)
        {
          ind = it->first;
	  i1 = offset+ind; i2 = offset+ind+udim; i3 = offset+ind+2*udim;
	  x_d = mesh->getCoord(it->second, 0).getVal();
	  y_d = mesh->getCoord(it->second, 1).getVal();
	  z_d = mesh->getCoord(it->second, 2).getVal();
	  norm = n[ind].norm2();
	  n[ind](0) = n[ind](0)/norm; n[ind](1) = n[ind](1)/norm; n[ind](2) = n[ind](2)/norm;
	  calculateTangentsFromNormal(n[ind], t1, t2);
	  if (Uinfo[c]->normalbc.find(b->second) != Uinfo[c]->normalbc.end() && Uinfo[c]->tangentbc.find(b->second) != Uinfo[c]->tangentbc.end())
	  {
	    r1 = t1(0)*q(i1) + t1(1)*q(i2) + t1(2)*q(i3) - Uinfo[c]->tangentbc[b->second]->Eval();
	    r2 = t2(0)*q(i1) + t2(1)*q(i2) + t2(2)*q(i3) - Uinfo[c]->tangentbc[b->second]->Eval();
	    r3 = q(i1)*n[ind](0) + q(i2)*n[ind](1) + q(i3)*n[ind](2)-Uinfo[c]->normalbc[b->second]->Eval();
	    setResidual(r1, i1);
	    setResidual(r2, i2);
	    setResidual(r3, i3);
	  }
	  else if (Uinfo[c]->normalbc.find(b->second) != Uinfo[c]->normalbc.end())
	  {
	    d1.clear(); d2.clear();
	    for (t_drq::iterator it=drq[i1].begin();it!=drq[i1].end();it++)
	    {
	      d1[it->first] = t1(0).getVal()*it->second;
	      d2[it->first] = t2(0).getVal()*it->second;
	    }
	    for (t_drq::iterator it=drq[i2].begin();it!=drq[i2].end();it++)
	    {
	      d1[it->first] += t1(1).getVal()*it->second;
	      d2[it->first] += t2(1).getVal()*it->second;
	    }
	    for (t_drq::iterator it=drq[i3].begin();it!=drq[i3].end();it++)
	    {
	      d1[it->first] += t1(2).getVal()*it->second;
	      d2[it->first] += t2(2).getVal()*it->second;
	    }
	    drq[i1].clear();
	    drq[i1] = d1;
	    drq[i2].clear();
	    drq[i2] = d2;
	    r1 = t1(0)*r[i1] + t1(1)*r[i2] + t1(2)*r[i3];
	    r2 = t2(0)*r[i1] + t2(1)*r[i2] + t2(2)*r[i3];
	    r3 = q(i1)*n[ind](0) + q(i2)*n[ind](1) + q(i3)*n[ind](2)
	         - Uinfo[c]->normalbc[b->second]->Eval();
	    r[i1] = r1.getVal();
	    r[i2] = r2.getVal();
	    setResidual(r3, i3);
	  }
	  else if (Uinfo[c]->tangentbc.find(b->second) != Uinfo[c]->tangentbc.end())
	  {
	    d1.clear();
	    for (t_drq::iterator it=drq[offset+ind].begin();it!=drq[offset+ind].end();it++) d1[it->first] = n[ind](0).getVal()*it->second;
	    for (t_drq::iterator it=drq[offset+ind+udim].begin();it!=drq[offset+ind+udim].end();it++) d1[it->first] += n[ind](1).getVal()*it->second;
	    for (t_drq::iterator it=drq[offset+ind+2*udim].begin();it!=drq[offset+ind+2*udim].end();it++) d1[it->first] += n[ind](2).getVal()*it->second;
	    drq[ind+offset+2*udim].clear();
	    drq[ind+offset+2*udim] = d1;
	    r1 = t1(0)*q(offset+ind) + t1(1)*q(offset+ind+udim) + t1(2)*q(offset+ind+2*udim) - Uinfo[c]->tangentbc[b->second]->Eval();
	    r2 = t2(0)*q(offset+ind) + t2(1)*q(offset+ind+udim) + t2(2)*q(offset+ind+2*udim) - Uinfo[c]->tangentbc[b->second]->Eval();
	    r3 = n[ind](0)*r[offset+ind] + n[ind](1)*r[offset+udim+ind] + n[ind](2)*r[offset+2*udim+ind];
	    r[ind+offset+2*udim] = r3.getVal();
	    setResidual(r1, ind+offset);
	    setResidual(r2, ind+offset+udim);
	  }
        }
      }
      offset += Uinfo[c]->type*udim;
    }
  }
}

void CStateProblem::assembleEdgeNormalTangentBC()
{
  int offset, c, e, i, j, ind;
  IntPoint<1> *ip;
  
  if (dim != 3) return;
  
  for (map<int,string>::iterator b=mesh->beginPhysicalNames(); b != mesh->endPhysicalNames(); b++)
  {
    offset = 0;
    for (c=0; c<(int)Uinfo.size(); c++)
    {
      ip = U[Uinfo[c]->index]->ip1d;
      int udim=ip->getMeshNdof();
      if (Uinfo[c]->normalbc.find(b->second) != Uinfo[c]->normalbc.end() || Uinfo[c]->tangentbc.find(b->second) != Uinfo[c]->tangentbc.end())
      {
        map<int,int> index;
        map<int,CVec> t;
	CVec n1, n2, lt;
        index.clear();
        t.clear();
        for (e=0; e<mesh->getNedges(); e++)
	{
	  if (mesh->getPhysicalName(mesh->getEdgePhysicalNo(e)) == b->second)
	  {
	    mesh->calculateEdgeTangent(e, lt);
	    for (i=0; i<ip->getNdof(); i++)
	    {
	      ind = ip->mapDof(e,i);
	      if (index.find(ind) == index.end())
	      {
	        index[ind] = mesh->getEdgeNode(e, i);
		t    [ind] = lt;
	      }
	      else
	      {
		for (j=0; j<3; j++) t[ind](j) += lt(j);
	      }
	    }
	  }
	}
	CVar r1, r2, r3;
	double norm;
	t_drq d1, d2;
        for (map<int,int>::iterator it=index.begin(); it != index.end(); it++)
        {
          ind = it->first;
	  x_d = mesh->getCoord(it->second, 0).getVal();
	  y_d = mesh->getCoord(it->second, 1).getVal();
	  z_d = mesh->getCoord(it->second, 2).getVal();
	  norm = t[ind].norm2();
	  t[ind](0) = t[ind](0)/norm; t[ind](1) = t[ind](1)/norm; t[ind](2) = t[ind](2)/norm;
	  calculateTangentsFromNormal(t[ind], n1, n2);
	  if (Uinfo[c]->normalbc.find(b->second) != Uinfo[c]->normalbc.end() && Uinfo[c]->tangentbc.find(b->second) != Uinfo[c]->tangentbc.end())
	  {
	    r1 = n1(0)*q(offset+ind) + n1(1)*q(offset+ind+udim) + n1(2)*q(offset+ind+2*udim)-Uinfo[c]->normalbc[b->second]->Eval();
	    r2 = n2(0)*q(offset+ind) + n2(1)*q(offset+ind+udim) + n2(2)*q(offset+ind+2*udim)-Uinfo[c]->normalbc[b->second]->Eval();
	    r3 = q(offset+ind)*t[ind](0) + q(offset+udim+ind)*t[ind](1) + q(offset+2*udim+ind)*t[ind](2)-Uinfo[c]->tangentbc[b->second]->Eval();
	    setResidual(r1, ind+offset);
	    setResidual(r2, ind+offset+udim);
	    setResidual(r3, ind+offset+2*udim);
	  }
	  else if (Uinfo[c]->tangentbc.find(b->second) != Uinfo[c]->tangentbc.end())
	  {
	    d1.clear(); d2.clear();
	    for (t_drq::iterator it=drq[offset+ind].begin();it!=drq[offset+ind].end();it++)
	    {
	      d1[it->first] = n1(0).getVal()*it->second;
	      d2[it->first] = n2(0).getVal()*it->second;
	    }
	    for (t_drq::iterator it=drq[offset+ind+udim].begin();it!=drq[offset+ind+udim].end();it++)
	    {
	      d1[it->first] += n1(1).getVal()*it->second;
	      d2[it->first] += n2(1).getVal()*it->second;
	    }
	    for (t_drq::iterator it=drq[offset+ind+2*udim].begin();it!=drq[offset+ind+2*udim].end();it++)
	    {
	      d1[it->first] += n1(2).getVal()*it->second;
	      d2[it->first] += n2(2).getVal()*it->second;
	    }
	    drq[ind+offset].clear();
	    drq[ind+offset] = d1;
	    drq[ind+offset+udim].clear();
	    drq[ind+offset+udim] = d2;
	    r1 = n1(0)*r[offset+ind] + n1(1)*r[offset+ind+udim] + n1(2)*r[offset+ind+2*udim];
	    r2 = n2(0)*r[offset+ind] + n2(1)*r[offset+ind+udim] + n2(2)*r[offset+ind+2*udim];
	    r3 = q(offset+ind)*t[ind](0) + q(offset+udim+ind)*t[ind](1) + q(offset+2*udim+ind)*t[ind](2)
	         - Uinfo[c]->tangentbc[b->second]->Eval();
	    r[ind+offset] = r1.getVal();
	    r[ind+offset+udim] = r2.getVal();
	    setResidual(r3, ind+offset+2*udim);
	  }
	  else if (Uinfo[c]->normalbc.find(b->second) != Uinfo[c]->normalbc.end())
	  {
	    d1.clear();
	    for (t_drq::iterator it=drq[offset+ind].begin();it!=drq[offset+ind].end();it++) d1[it->first] = t[ind](0).getVal()*it->second;
	    for (t_drq::iterator it=drq[offset+ind+udim].begin();it!=drq[offset+ind+udim].end();it++) d1[it->first] += t[ind](1).getVal()*it->second;
	    for (t_drq::iterator it=drq[offset+ind+2*udim].begin();it!=drq[offset+ind+2*udim].end();it++) d1[it->first] += t[ind](2).getVal()*it->second;
	    drq[ind+offset+2*udim].clear();
	    drq[ind+offset+2*udim] = d1;
	    r1 = n1(0)*q(offset+ind) + n1(1)*q(offset+ind+udim) + n1(2)*q(offset+ind+2*udim) - Uinfo[c]->normalbc[b->second]->Eval();
	    r2 = n2(0)*q(offset+ind) + n2(1)*q(offset+ind+udim) + n2(2)*q(offset+ind+2*udim) - Uinfo[c]->normalbc[b->second]->Eval();
	    r3 = t[ind](0)*r[offset+ind] + t[ind](1)*r[offset+udim+ind] + t[ind](2)*r[offset+2*udim+ind];
	    r[ind+offset+2*udim] = r3.getVal();
	    setResidual(r1, ind+offset);
	    setResidual(r2, ind+offset+udim);
	  }
        }
      }
      offset += Uinfo[c]->type*udim;
    }
  }
}

void CStateProblem::setDirichletBC()
{
	int c, e, i, ind, offset;
	log(W_SOLVE, "-Setting boundary conditions.\n");
  
	for (e=0; e<mesh->getNvolumes(); e++)
	{
		IntPoint<3> *ip;
		Type<3>::t_IP::iterator it;
		for (it=IP3D.begin(); it!=IP3D.end(); it++)
		{
		  ip   = it->second;
		  for (i=0; i<min(ip->getNdof(),mesh->getNvolumeDofs()); i++)
		  {
			ind = ip->mapDof(e, i);
			x = mesh->getCoord(mesh->getVolumeNode(e,i),0);
			y = mesh->getCoord(mesh->getVolumeNode(e,i),1);
			z = mesh->getCoord(mesh->getVolumeNode(e,i),2);
			x_d = x.getVal();
			y_d = y.getVal();
			z_d = z.getVal();
			offset = 0;
			for (c=0; c<(int)U.size(); c++)
			{
			  if (!U[c]->ip3d->isDiscontinuous() && U[c]->ip3d == ip && U[c]->p.find(mesh->getPhysicalName(mesh->getVolumePhysicalNo(e))) != U[c]->p.end())
			  {
				q(offset+ind) = U[c]->p[mesh->getPhysicalName(mesh->getVolumePhysicalNo(e))]->Eval();
			  }
			  offset += U[c]->ip3d->getMeshNdof();
			}
		  }
		}
	}
  
	for (e=0; e<mesh->getNfaces(); e++)
	{
		IntPoint<2> *ip;
		Type<2>::t_IP::iterator it;
		for (it=IP2D.begin(); it!=IP2D.end(); it++)
		{
		  ip   = it->second;
		  for (i=0; i<min(ip->getNdof(),mesh->getNfaceDofs()); i++)
		  {
			ind = ip->mapDof(e, i);
			x = mesh->getCoord(mesh->getFaceNode(e,i),0);
			y = mesh->getCoord(mesh->getFaceNode(e,i),1);
			z = mesh->getCoord(mesh->getFaceNode(e,i),2);
			x_d = x.getVal();
			y_d = y.getVal();
			z_d = z.getVal();
			offset = 0;
			for (c=0; c<(int)U.size(); c++)
			{
			  if (!U[c]->ip2d->isDiscontinuous() && U[c]->ip2d == ip && U[c]->p.find(mesh->getPhysicalName(mesh->getFacePhysicalNo(e))) != U[c]->p.end())
			  {
				q(offset+ind) = U[c]->p[mesh->getPhysicalName(mesh->getFacePhysicalNo(e))]->Eval();
			  }
			  offset += U[c]->ip2d->getMeshNdof();
			}
		  }
		}
	}
	for (e=0; e<mesh->getNedges(); e++)
	{
		IntPoint<1> *ip;
		Type<1>::t_IP::iterator it;
		for (it=IP1D.begin(); it!=IP1D.end(); it++)
		{
		  ip   = it->second;
		  for (i=0; i<min(ip->getNdof(),mesh->getNedgeDofs()); i++)
		  {
			ind = ip->mapDof(e, i);
			x = mesh->getCoord(mesh->getEdgeNode(e,i),0);
			y = mesh->getCoord(mesh->getEdgeNode(e,i),1);
			z = mesh->getCoord(mesh->getEdgeNode(e,i),2);
			x_d = x.getVal();
			y_d = y.getVal();
			z_d = z.getVal();
			offset = 0;
			for (c=0; c<(int)U.size(); c++)
			{
			  if (!U[c]->ip1d->isDiscontinuous() && U[c]->ip1d == ip && U[c]->p.find(mesh->getPhysicalName(mesh->getEdgePhysicalNo(e))) != U[c]->p.end())
			  {
				q(offset+ind) = U[c]->p[mesh->getPhysicalName(mesh->getEdgePhysicalNo(e))]->Eval();
			  }
			  offset += U[c]->ip1d->getMeshNdof();
			}
		  }
		}
	}
	for (e=0; e<mesh->getNpoints(); e++)
	{
		ind = mesh->getPointNode(e);
		x = mesh->getCoord(ind,0);
		y = mesh->getCoord(ind,1);
		z = mesh->getCoord(ind,2);
		x_d = x.getVal();
		y_d = y.getVal();
		z_d = z.getVal();
		offset = 0;
		for (c=0; c<(int)U.size(); c++)
		{
		  if (!U[c]->ip1d->isDiscontinuous() && U[c]->p.find(mesh->getPhysicalName(mesh->getPointPhysicalNo(e))) != U[c]->p.end())
		  {
			q(offset+ind) = U[c]->p[mesh->getPhysicalName(mesh->getPointPhysicalNo(e))]->Eval();
		  }
		  offset += U[c]->ip1d->getMeshNdof();
		}
	}
}

void CStateProblem::setResidualsForDirichletBC()
{
	int c, e, i, ind, offset;
  
	for (e=0; e<mesh->getNfaces(); e++)
	{
		IntPoint<2> *ip;
		Type<2>::t_IP::iterator it;
		for (it=IP2D.begin(); it!=IP2D.end(); it++)
		{
		  ip   = it->second;
		  for (i=0; i<ip->getNdof(); i++)
		  {
			ind = ip->mapDof(e, i);
			offset = 0;
			for (c=0; c<(int)U.size(); c++)
			{
			  if (!U[c]->ip2d->isDiscontinuous() && U[c]->ip2d == ip && U[c]->p.find(mesh->getPhysicalName(mesh->getFacePhysicalNo(e))) != U[c]->p.end())
			  {
				r[offset+ind] = 0;
				drq[offset+ind].clear();
				drq[offset+ind][offset+ind] = 1;
			  }
			  offset += U[c]->ip2d->getMeshNdof();
			}
		  }
		}
  	}
	for (e=0; e<mesh->getNedges(); e++)
	{
		IntPoint<1> *ip;
		Type<1>::t_IP::iterator it;
		for (it=IP1D.begin(); it!=IP1D.end(); it++)
		{
		  ip   = it->second;
		  for (i=0; i<ip->getNdof(); i++)
		  {
			ind = ip->mapDof(e, i);
			offset = 0;
			for (c=0; c<(int)U.size(); c++)
			{
			  if (!U[c]->ip1d->isDiscontinuous() && U[c]->ip1d == ip && U[c]->p.find(mesh->getPhysicalName(mesh->getEdgePhysicalNo(e))) != U[c]->p.end())
			  {
					r[offset+ind] = 0;
					drq[offset+ind].clear();
					drq[offset+ind][offset+ind] = 1;
			  }
			  offset += U[c]->ip1d->getMeshNdof();
			}
		  }
		}
	}
	for (e=0; e<mesh->getNpoints(); e++)
	{
		ind = mesh->getPointNode(e);
		offset = 0;
		for (c=0; c<(int)U.size(); c++)
		{
		  if (!U[c]->ip1d->isDiscontinuous() && U[c]->p.find(mesh->getPhysicalName(mesh->getPointPhysicalNo(e))) != U[c]->p.end())
		  {
			r[offset+ind] = 0;
			drq[offset+ind].clear();
			drq[offset+ind][offset+ind] = 1;
		  }
		  offset += U[c]->ip1d->getMeshNdof();
		}
	}
}

void CStateProblem::setMeanValues()
{
  int c, i, e, offset = 0, ind0, ind;
  for (c=0; c<(int)U.size(); c++)
  {
    if (U[c]->set_mean_value)
    {
      ind0 = -1;
      for (e=0; e<mesh->getNvolumes(); e++)
      {
        if (U[c]->mv_segment == mesh->getVolumePhysicalNo(e))
	{
	  for (i=0; i<U[c]->ip3d->getNdof(); i++)
	  {
	    ind = U[c]->ip3d->mapDof(e, i);
	    if (ind0 == -1)
	    {
              ind0 = ind;
              r[offset+ind0] = -U[c]->mean_value*U[c]->volume;
              drq[offset+ind0].clear();
	    }
	    if (U[c]->mean[ind] != 0)
	    {
              r[offset+ind0] += q(offset+ind).getVal()*U[c]->mean[ind];
              drq[offset+ind0][offset+ind] = U[c]->mean[ind];
	      U[c]->mean[ind] = 0;
	    }
	  }
	}
      }
      for (e=0; e<mesh->getNfaces(); e++)
      {
        if (U[c]->mv_segment == mesh->getFacePhysicalNo(e))
	{
	  for (i=0; i<U[c]->ip2d->getNdof(); i++)
	  {
	    ind = U[c]->ip2d->mapDof(e, i);
	    if (ind0 == -1)
	    {
	      ind0 = ind;
              r[offset+ind0] = -U[c]->mean_value*U[c]->volume;
              drq[offset+ind0].clear();
	    }
	    if (U[c]->mean[ind] != 0)
	    {
              r[offset+ind0] += q(offset+ind).getVal()*U[c]->mean[ind];
              drq[offset+ind0][offset+ind] = U[c]->mean[ind];
	      U[c]->mean[ind] = 0;
	    }
	  }
	}
      }
      for (e=0; e<mesh->getNedges(); e++)
      {
        if (U[c]->mv_segment == mesh->getEdgePhysicalNo(e))
	{
	  for (i=0; i<U[c]->ip1d->getNdof(); i++)
	  {
	    ind = U[c]->ip1d->mapDof(e, i);
	    if (ind0 == -1)
	    {
	      ind0 = ind;
              r[offset+ind0] = -U[c]->mean_value*U[c]->volume;
              drq[offset+ind0].clear();
	    }
	    if (U[c]->mean[ind] != 0)
	    {
              r[offset+ind0] += q(offset+ind).getVal()*U[c]->mean[ind];
              drq[offset+ind0][offset+ind] = U[c]->mean[ind];
	      U[c]->mean[ind] = 0;
	    }
	  }
	}
      }
    }
    offset += U[c]->ip3d->getMeshNdof();
  }
}

void CStateProblem::assemble()
{
  clock_t t;
  
  log(W_SOLVE, "Starting assembly:\n");
  start_clock(t);
  
  for (int i=0; i<q.size(); i++)
  {
    r[i] = 0;
    drq[i].clear();
  }
  for (int c=0; c<(int)U.size(); c++)
  {
    if (!U[c]->set_mean_value || U[c]->mean == 0) continue;
    for (int i=0; i<U[c]->ip3d->getMeshNdof(); i++)
    {
      U[c]->mean[i] = 0;
    }
    U[c]->volume = 0;
  }
  
  // set boundary conditions
  setDirichletBC();

  // assemble integrals
  if (useVolumeIntegral()) assembleVolumeIntegral();
  if (useFaceIntegral()) assembleFaceIntegral();
  if (useEdgeIntegral()) assembleEdgeIntegral();
  
  // set normal and tangent boundary conditions
  assembleFaceNormalTangentBC();
  assembleEdgeNormalTangentBC();

  // set boundary conditions
  setResidualsForDirichletBC();

  // set mean values
  setMeanValues();
  
  afterAssemble();

  log(W_SOLVE, "Assembly completed in %f s.\n", stop_clock(t));
}

void CStateProblem::createSparseMatAndRhs(MatrixUMF &sp, double *&f, vector<int> &indices)
{
	  t_drq::iterator it;
	  int i, j, size=0;
	  int nnz = 0;
	  int *ia, *ja;
	  double *v;
	  vector<int> inv_indices;

	  inv_indices.resize(q.size(), -1);
	  for (i=0; i<q.size(); i++)
	  {
	    if (drq[i].size() > 0)
	    {
	      inv_indices[i] = size;
	      size++;
	      indices.push_back(i);
	      nnz += drq[i].size();
	    }
	  }
	  sp.resize(size, size, nnz);
	  ia = sp.getIa();
	  ja = sp.getJa();
	  v  = sp.getV();
	  f = new double[size];
	  for (i=0; i<size; i++)
	  {
	    ia[i+1] = ia[i] + drq[indices[i]].size();
	    for (j=0, it=drq[indices[i]].begin(); it !=drq[indices[i]].end(); j++, it++)
	    {
	        ja[ia[i]+j] = inv_indices[it->first];
	        v[ia[i]+j]  = it->second;
	    }
	    f[i] = r[indices[i]];
	    drq[indices[i]].clear();
	  }
}

void CStateProblem::createPetscMatAndRhs(Mat &M, Vec &f)
{
	int *cols, *rows, nnz, nrows, offset;
	double *vals;
	bool active_rows[q.size()]; // true for rows belonging to the actual partition.

	// count nonzero elements in system matrix
	nnz = 0;
	nrows = 0;
	offset = 0;
	for (unsigned int i=0; i<U.size(); i++)
	{
		for (int j=0; j<U[i]->ip3d->getMeshNdof(); j++)
		{
			if (dof_distribution->np() == 1 || mesh->getNodePartNo(j) == dof_distribution->myp()+1)
			{
				nnz += drq[offset+j].size();
				nrows++;
				active_rows[offset+j] = true;
			}
			else
			{
				drq[offset+j].clear();
				active_rows[offset+j] = false;
			}
		}
		offset += U[i]->ip3d->getMeshNdof();
	}

	ASSERT(nrows == dof_distribution->lsize(), "Number of rows mismatch: %d != %d!", nrows, dof_distribution->lsize());

	vals = new double[nnz];
	cols = new int[nnz];
	rows = new int[nrows+1];

	// create CSR arrays
	int n = 0, irow = 0;
	int p = (dof_distribution->np() == 1)?0:dof_distribution->myp()+1;
	offset = 0;
	for (unsigned int c=0; c<U.size(); c++)
	{
		for (set<int>::iterator i=mesh->getPartitionNodes(p).begin(); i!=mesh->getPartitionNodes(p).end(); i++)
		{
			rows[irow] = n;
			for (t_drq::iterator j=drq[offset+*i].begin(); j!=drq[offset+*i].end(); j++)
			{
				vals[n] = j->second;
				cols[n] = petsc_dof_numbering[j->first];
				n++;
			}
			drq[offset+*i].clear();
			irow++;
		}
		offset += U[c]->ip3d->getMeshNdof();
	}

	rows[nrows] = n;

	// create MPI matrix
	MatCreateMPIAIJWithArrays(PETSC_COMM_WORLD, dof_distribution->lsize(), dof_distribution->lsize(), PETSC_DETERMINE, PETSC_DETERMINE, rows, cols, vals, &M);

	delete[] vals;
	delete[] cols;
	delete[] rows;

	// create vector of rhs
	int j=0;
	for (int i=0; i<q.size(); i++)
	{
		if (active_rows[i])
		{
			f_array[j] = r[i];
			j++;
		}
	}
	VecCreateMPIWithArray(PETSC_COMM_WORLD, dof_distribution->lsize(), PETSC_DECIDE, f_array, &f);
}

int CStateProblem::solve(double &xnorm, double relax)
{
  MatrixUMF sp;
  int status;
  double *f = 0;
  clock_t t;
  vector<int> indices; //< Indices of nonzero matrix rows.

  log(W_SOLVE, "Allocating sparse matrix.\n");
  createSparseMatAndRhs(sp, f, indices);

  log(W_SOLVE, "Starting solution.\n");
  start_clock(t);
  status = sp.solve(f);

  if (status < 0)
  {
    log(W_SOLVE, "Solution failed (error %d).\n", status);
    return status;
  }

  xnorm = 0;
  for (unsigned int i=0; i<indices.size(); i++)
  {
    xnorm += f[i]*f[i];
    q(indices[i]) = q(indices[i])-f[i]*relax;
    r[indices[i]] = f[i];
  }
  xnorm = sqrt(xnorm);
  delete[] f;

  log(W_SOLVE, "Solution completed in %f s.\n\n", stop_clock(t));

  return 0;
}

int CStateProblem::solvePetsc(double &xnorm, double relax)
{
  Mat A;
  Vec b;
  KSP ksp;
  int status;
  clock_t t;

  log(W_SOLVE, "Allocating sparse matrix.\n");
  createPetscMatAndRhs(A, b);
  
  log(W_SOLVE, "Starting solution.\n");
  start_clock(t);
  KSPCreate(PETSC_COMM_WORLD, &ksp);
  KSPSetOperators(ksp, A, A, DIFFERENT_NONZERO_PATTERN);
  KSPSetFromOptions(ksp);
  status = KSPSolve(ksp, b, b);
  
  if (status != 0)
  {
    log(W_SOLVE, "Solution failed (error %d).\n", status);
    return status;
  }
  
  xnorm = 0;
  int offset = 0, ind, k = 0;
  for (unsigned int i=0; i<U.size(); i++)
  {
	  for (int j=0; j<U[i]->ip3d->getMeshNdof(); j++)
	  {
		  if (dof_distribution->np() == 1 || mesh->getNodePartNo(j) == dof_distribution->myp()+1)
		  {
			  ind = offset + j;
			  q(ind) = q(ind)-f_array[k]*relax;
			  k++;
		  }
	  }
	  offset += U[i]->ip3d->getMeshNdof();
  }
  VecNorm(b, NORM_2, &xnorm);
  KSPDestroy(&ksp);
  VecDestroy(&b);
  MatDestroy(&A);

  // update solution on nodes shared by multiple partitions
  const int sd_size = U.size()*mesh->getSharedNodes().size();
  double shared_dofs_out[sd_size], shared_dofs_in[dof_distribution->np()*sd_size];
  offset = 0;
  int j=0;
  for (unsigned int c=0; c<U.size(); c++)
  {
	  for (set<int>::iterator i=mesh->getSharedNodes().begin(); i!=mesh->getSharedNodes().end(); i++)
	  {
		  shared_dofs_out[j] = q(offset+*i).getVal();
		  j++;
	  }
	  offset += U[c]->ip3d->getMeshNdof();
  }
  MPI_Allgather(shared_dofs_out, sd_size, MPI_DOUBLE, shared_dofs_in, sd_size, MPI_DOUBLE, PETSC_COMM_WORLD);
  offset = 0;
  int poffset = (dof_distribution->np() == 1)?0:1;
  j=0;
  for (unsigned int c=0; c<U.size(); c++)
  {
	  for (set<int>::iterator i=mesh->getSharedNodes().begin(); i!=mesh->getSharedNodes().end(); i++)
	  {
		  if (mesh->getNodePartNo(*i) != dof_distribution->myp()+poffset)
			  q(offset+*i) = q(offset+*i) - q(offset+*i).getVal() + shared_dofs_in[(mesh->getNodePartNo(*i)-poffset)*sd_size + j];
		  j++;
	  }
	  offset += U[c]->ip3d->getMeshNdof();
  }



  log(W_SOLVE, "Solution completed in %f s.\n\n", stop_clock(t));
    
  return 0;
}

void CStateProblem::nlSolve(double relax)
{
  int it = 0, status = 0;
  clock_t t;
  double deltanorm, t_asm, t_sol;

  log(W_PROG, "\nStarting nonlinear loop.\n");
  log(W_PROG, "+----+---------+---------+--------+--------+\n");
  log(W_PROG, "|iter|residuum |deltanorm|time_asm|time_sol|\n");
  log(W_PROG, "+----+---------+---------+--------+--------+\n");
  
  start_clock(t);
  
  deltanorm = 1e30;
  rnorm = 1e30;
  
  assemble();
  t_asm = stop_clock(t);
  
  // calculate residual norm
  double lrnorm = 0;
  int p = (dof_distribution->np() == 1)?0:dof_distribution->myp()+1;
  int offset = 0;
  for (unsigned int c=0; c<U.size(); c++)
  {
	  for (set<int>::iterator i=mesh->getPartitionNodes(p).begin(); i!=mesh->getPartitionNodes(p).end(); i++)
		  lrnorm += r[offset+*i]*r[offset+*i];
	  offset += U[c]->ip3d->getMeshNdof();
  }
  lrnorm = sqrt(lrnorm);
  MPI_Allreduce(&lrnorm, &rnorm, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
  
  try
  {
    while (rnorm > rstop)
    {
      start_clock(t);
      status = solvePetsc(deltanorm, 1.);
      if (status < 0) throw status;
      afterNonlinearStep(it);
    
      t_sol = stop_clock(t);
    
      log(W_PROG, "|%4d|%-9.3e|%-9.3e|%8.2f|%8.2f|\n", it, rnorm, deltanorm, t_asm, t_sol);
    
      it++;
      log(W_SOLVE, "Nonlinear iteration no. %d:\n", it);
    
      start_clock(t);
      assemble();
      t_asm = stop_clock(t);
    
      // calculate residual norm
      lrnorm = 0;
      offset = 0;
      for (unsigned int c=0; c<U.size(); c++)
      {
    	  for (set<int>::iterator i=mesh->getPartitionNodes(p).begin(); i!=mesh->getPartitionNodes(p).end(); i++)
    		  lrnorm += r[offset+*i]*r[offset+*i];
    	  offset += U[c]->ip3d->getMeshNdof();
      }
      lrnorm = sqrt(lrnorm);
      MPI_Allreduce(&lrnorm, &rnorm, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
    }
  }
  catch (int e)
  {
    log(W_PROG, "|%4d|Solution error (%2d)|%8.2f|  ----  |\n", it, status, t_asm);
    log(W_PROG, "+----+---------+---------+--------+--------+\n\n");
    return;
  }
  
  log(W_PROG, "|%4d|%-9.3e|  -----  |%8.2f|  ----  |\n", it, rnorm, t_asm);
  log(W_PROG, "+----+---------+---------+--------+--------+\n\n");
}

void CStateProblem::tSolve(double relax)
{
  int it = 0;
  clock_t t;
  double deltanorm, t_asm, t_sol;
  char fn[255], pvd_fn[255];

  log(W_PROG, "\nStarting time loop.\n");
  log(W_PROG, "+--------+--------+--------+\n");
  log(W_PROG, "|time_sim|time_asm|time_sol|\n");
  log(W_PROG, "+--------+--------+--------+\n");
  
  oq.resize(q.size());
  ooq.resize(q.size());
  sprintf(fn, "%s.%d.vtu", params["savebase"].c_str(), it);
  
  int pos = params["savebase"].find_last_of('/');
  if (pos > (int)params["savebase"].length()) pos = 0;
  if (dof_distribution->np() == 1)
  {
	  sprintf(pvd_fn, "%s.pvd", params["savebase"].substr(pos+1).c_str());
  }
  else
  {
	  sprintf(pvd_fn, "%s.np%d.pvd", params["savebase"].substr(pos+1).c_str(), dof_distribution->np());
  }
  ofstream f;

  //saveVTKview(fn);
  if (dof_distribution->myp() == 0)
  {
	  f.open(pvd_fn);
	  f << "<?xml version=\"1.0\"?>\n"
	    << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n"
        << "<Collection>\n";
//        << "<DataSet timestep=\"0\" group=\"\" part=\"0\" file=\"" << params["savebase"].c_str() << ".0.vtu\"/>\n";
  }
  
  while (tau*it < tmax)
  {
    it++;
    time = tau*it;
    log(W_SOLVE, "Solving time %.2f:\n", time);

    for (int i=0; i<q.size(); ooq(i) = oq(i).getVal(), oq(i) = q(i).getVal(), i++);
    calculated_solution = false;
    
    start_clock(t);
    assemble();
    t_asm = stop_clock(t);
    
    // calculate residual norm
    double lrnorm = 0;
    int p = (dof_distribution->np() == 1)?0:dof_distribution->myp()+1;
    int offset = 0;
    for (unsigned int c=0; c<U.size(); c++)
    {
  	  for (set<int>::iterator i=mesh->getPartitionNodes(p).begin(); i!=mesh->getPartitionNodes(p).end(); i++)
  		  lrnorm += r[offset+*i]*r[offset+*i];
  	  offset += U[c]->ip3d->getMeshNdof();
    }
    lrnorm = sqrt(lrnorm);
    MPI_Allreduce(&lrnorm, &rnorm, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);

    start_clock(t);
    solvePetsc(deltanorm, relax);
    t_sol = stop_clock(t);
    
    log(W_PROG, "|%8.2f|%8.2f|%8.2f|\n", time, t_asm, t_sol);
    
    if (dof_distribution->np() == 1)
    {
    	sprintf(fn, "%s.%d.vtu", params["savebase"].c_str(), it);
    }
    else
    {
    	sprintf(fn, "%s.%dof%d.%d.vtu", params["savebase"].c_str(), dof_distribution->myp()+1, dof_distribution->np(), it);
    }
    saveVTKview(fn);
    if (dof_distribution->myp() == 0)
    {
    	for (int i=1; i<=dof_distribution->np(); i++)
    	{
    		if (dof_distribution->np() == 1)
			{
				sprintf(fn, "%s.%d.vtu", params["savebase"].c_str(), it);
			}
			else
			{
				sprintf(fn, "%s.%dof%d.%d.vtu", params["savebase"].c_str(), i, dof_distribution->np(), it);
			}
    		f << "<DataSet timestep=\"" << it << "\" group=\"\" part=\"" << i << "\" file=\"" << fn << "\"/>\n";
    	}
    	f.flush();
    }

    sprintf(fn, "%s.%dof%d.%d.sol", params["savebase"].c_str(), dof_distribution->myp()+1, dof_distribution->np(), it);
    saveSolution(fn);
  }
  
  log(W_PROG, "+--------+--------+--------+\n\n");
  
  if (dof_distribution->myp() == 0)
  {
	  f << "</Collection>\n"
        << "</VTKFile>";
	  f.close();
  }
  log(W_SOLVE, "Animation saved to %s.", fn);
}

void CStateProblem::saveMatrix(const char *file)
{
  int n = q.size();
  ofstream f(file);

  for (int i=0; i<n; i++)
  {
    for (t_drq::iterator it=drq[i].begin(); it!=drq[i].end(); it++)
    {
      if (it->second != 0)
        f << i << " " << it->first  << " " << it->second << endl;
    }
  }
  
  f.close();
}

void CStateProblem::addPPscalar(string name, string expression)
{
  vector<string> expr;
  expr.push_back(expression);
  PPfields[name] = new PPdata(name, UT_SCALAR, expr, mesh->getNvolumes()*mesh->getNvolumeDofs()+mesh->getNfaces()*mesh->getNfaceDofs()+mesh->getNedges()*mesh->getNedgeDofs());
  int d = output->registerNodeData(name, UT_SCALAR);
  output->addComponent(d, PPfields[name]->getData(0));
}


void CStateProblem::calculateSolution()
{
  int e, k, j, n, c, offset, ind;
  CVar val;
  CVec grad;
  
  if (calculated_solution) return;
  
  for (int i=0; i<q.size(); i++)
  {
    if (isnan(q(i).getVal()))
    {
      q(i) = 1e30;
    }
  }
  
  int N = output->size();
  
  if (output->saveCells(Save3D))
  {
	  for (e=0; e<mesh->getNvolumes(); e++)
	  {
	    	// skip volumes of different partitions
	    	if (dof_distribution->np() > 1 && mesh->getVolumePartNo(e) != dof_distribution->myp()+1) continue;

		for (Type<3>::t_IP::iterator it=IP3D.begin(); it != IP3D.end(); it++)
		{
		  IntPoint<3> *ip = it->second;
		  ip->calcE(e, EVAL_FG);
		  for (k=0; k<mesh->getNvolumeDofs(); k++)
		  {
			ip->calcQP(k+ip->getQ()->n, EVAL_FG);
			ind = e*mesh->getNvolumeDofs() + k;//mesh->getVolumeNode(e, k);
			offset = 0;
			for (c=0; c<(int)U.size(); c++)
			{
			  if (U[c]->ip3d == ip)
			  {
				ip->eval(q, e, offset, val, grad, EVAL_FG);
				U[c]->node_value[ind] = val.getVal();
				U[c]->node_grad[    ind] = grad(0).getVal();
				U[c]->node_grad[N  +ind] = grad(1).getVal();
				U[c]->node_grad[N*2+ind] = grad(2).getVal();
			  }
			  offset += U[c]->ip3d->getMeshNdof();
			}
		  }
		}
	  }
  }
  
  if (output->saveCells(Save2D))
  {
	  for (e=0; e<mesh->getNfaces(); e++)
	  {
	    	// skip volumes of different partitions
	    	if (dof_distribution->np() > 1 && mesh->getFacePartNo(e) != dof_distribution->myp()+1) continue;

		for (Type<2>::t_IP::iterator it=IP2D.begin(); it != IP2D.end(); it++)
		{
		  IntPoint<2> *ip = it->second;
		  ip->calcE(e, EVAL_FG);
		  for (k=0; k<mesh->getNfaceDofs(); k++)
		  {
			ip->calcQP(k+ip->getQ()->n, EVAL_FG);
			ind = output->saveCells(Save3D)*mesh->getNvolumes()*mesh->getNvolumeDofs() + e*mesh->getNfaceDofs() + k;//mesh->getFaceNode(e, k);
			offset = 0;
			for (c=0; c<(int)U.size(); c++)
			{
			  if (U[c]->ip2d == ip)
			  {
				ip->eval(q, e, offset, val, grad, EVAL_FG);
				U[c]->node_value[ind] = val.getVal();
				U[c]->node_grad[    ind] = grad(0).getVal();
				U[c]->node_grad[N  +ind] = grad(1).getVal();
				U[c]->node_grad[N*2+ind] = grad(2).getVal();
			  }
			  offset += U[c]->ip2d->getMeshNdof();
			}
		  }
		}
	  }
  }
  
  if (output->saveCells(Save1D))
  {
	  for (e=0; e<mesh->getNedges(); e++)
	  {
	    	// skip volumes of different partitions
	    	if (dof_distribution->np() > 1 && mesh->getEdgePartNo(e) != dof_distribution->myp()+1) continue;

		for (Type<1>::t_IP::iterator it=IP1D.begin(); it != IP1D.end(); it++)
		{
		  IntPoint<1> *ip = it->second;
		  ip->calcE(e, EVAL_FG);
		  for (k=0; k<mesh->getNedgeDofs(); k++)
		  {
			ip->calcQP(k+ip->getQ()->n, EVAL_FG);
			ind = output->saveCells(Save3D)*mesh->getNvolumes()*mesh->getNvolumeDofs()
				+ output->saveCells(Save2D)*mesh->getNfaces()*mesh->getNfaceDofs()
				+ e*mesh->getNedgeDofs() + k;//mesh->getFaceNode(e, k);
			offset = 0;
			for (c=0; c<(int)U.size(); c++)
			{
			  if (U[c]->ip1d == ip)
			  {
				ip->eval(q, e, offset, val, grad, EVAL_FG);
				U[c]->node_value[ind] = val.getVal();
				U[c]->node_grad[    ind] = grad(0).getVal();
				U[c]->node_grad[N  +ind] = grad(1).getVal();
				U[c]->node_grad[N*2+ind] = grad(2).getVal();
			  }
			  offset += U[c]->ip3d->getMeshNdof();
			}
		  }
		}
	  }
  }
  

  char v[255];
  double node;
  
  sprintf(v, "node");
  
  // save additional scalar fields
  for (int i=0; i<N; i++)
  {
    node = i;
    
    Parser *p = new Parser;
    p->DefineVar(v, &node);

    // define parser variables
// the following is not working - coorditanes should be calculated in calculateSolution()
/*      x_d = mesh->getCoord(i, 0).getVal();
      y_d = mesh->getCoord(i, 1).getVal();
      z_d = mesh->getCoord(i, 2).getVal();
      p->DefineVar("x", &x_d);
      p->DefineVar("y", &y_d);
      p->DefineVar("z", &z_d);
//*/
    for (vector<CUnknown *>::iterator ui = Uinfo.begin(); ui != Uinfo.end(); ui++)
    {
      switch ((*ui)->type)
      {
        case UT_SCALAR:
	    p->DefineVar((*ui)->name, &U[(*ui)->index]->node_value[i]);
	    for (j=0; j<dim; j++)
            {
	      sprintf(v, "g%s%d", (*ui)->name.c_str(), j);
	      p->DefineVar(v, &U[(*ui)->index]->node_grad[N*j+i]);
	    }
	    break;
	  case UT_VECTOR:
	    for (int c=0; c<dim; c++)
	    {
	      sprintf(v, "%s%d", (*ui)->name.c_str(), c);
	      p->DefineVar(v, &U[(*ui)->index+c]->node_value[i]);
	      for (int j=0; j<dim; j++)
	      {
	        sprintf(v, "g%s%d%d", (*ui)->name.c_str(), c, j);
	        p->DefineVar(v, &U[(*ui)->index+c]->node_grad[N*j+i]);
	      }
	    }
	    break;
	  case UT_SYM_TENSOR:
	    const int *irow, *icol;
	    if (dim == 2)
	    {
	      n = 3;
	      irow = st2d_irow;
	      icol = st2d_icol;
	    }
	    else
	    {
	      n = 6;
	      irow = st_irow;
	      icol = st_icol;
	    }
	    for (int ti=0; ti<n; ti++)
	    {
	      sprintf(v, "%s%d%d", (*ui)->name.c_str(), irow[ti], icol[ti]);
	      p->DefineVar(v, &U[(*ui)->index+ti]->node_value[i]);
	      for (int j=0; j<dim; j++)
	      {
	        sprintf(v, "g%s%d%d%d", (*ui)->name.c_str(), irow[ti], icol[ti], j);
	        p->DefineVar(v, &U[(*ui)->index+ti]->node_grad[N*j+i]);
	      }
	      if (irow[ti] != icol[ti])
	      {
	        sprintf(v, "%s%d%d", (*ui)->name.c_str(), icol[ti], irow[ti]);
	        p->DefineVar(v, &U[(*ui)->index+ti]->node_value[i]);
	        for (int j=0; j<dim; j++)
	        {
	          sprintf(v, "g%s%d%d%d", (*ui)->name.c_str(), icol[ti], irow[ti], j);
	          p->DefineVar(v, &U[(*ui)->index+ti]->node_grad[N*j+i]);
	        }
	      }
	    }
	    break;
	  case UT_TENSOR:
	    for (int tr=0; tr<dim; tr++)
	    {
	      for (int tc=0; tc<dim; tc++)
	      {
	        sprintf(v, "%s%d%d", (*ui)->name.c_str(), tr, tc);
	        p->DefineVar(v, &U[(*ui)->index+tr*dim+tc]->node_value[i]);
	        for (int j=0; j<dim; j++)
	        {
	          sprintf(v, "g%s%d%d%d", (*ui)->name.c_str(), tr, tc, j);
	          p->DefineVar(v, &U[(*ui)->index+tr*dim+tc]->node_grad[N*j+i]);
	        }
	      }
	    }
	    break;
      }
    }

    for (map<string,PPdata*>::iterator it = PPfields.begin(); it != PPfields.end(); it++)
    {
      for (int c=0; c<it->second->getType(); c++)
      {
        p->SetExpr(it->second->getExpr(c));
        it->second->getData(c)[i] = p->Eval();
      }
    }

    delete p;
  }
  
  
  calculated_solution = true;
}


void CStateProblem::saveVTKview(const char *file)
{
  char fn[255], pn[255];

  if (file == 0)
  {
    if (params.find("savebase") != params.end())
    {
    	  if (dof_distribution->np() == 1)
    	  {
    		  sprintf(pn, "%s", "");
    	  }
    	  else
    	  {
    		  sprintf(pn,  ".p%dof%d", (dof_distribution->myp()+1), dof_distribution->np());
    	  }

    	  sprintf(fn, "%s%s.vtu", params["savebase"].c_str(), pn);
    	  file = fn;
    }
    else
    {
      log(tau>0?W_SOLVE:W_PROG, "VTK XML view will not be saved.\n");
      return;
    }
  }
  
  calculateSolution();
  
  output->saveVTK(file);
  
  log(tau>0?W_SOLVE:W_PROG, "Results saved to VTK file %s.\n", file);
}

bool CStateProblem::loadSolution(const char *fname)
{
	char fn[255], pn[255];
	if (fname == 0)
	{
	  if (params.find("loadfile") != params.end())
	  {
    	  fname = params["loadfile"].c_str();
	  }
	  else
	  {
		  loaded_solution = false;
		  return false;
	  }
	}

	if (dof_distribution->np() == 1)
	{
	  sprintf(pn, "%s", "");
	}
	else
	{
	  sprintf(pn,  ".p%dof%d", (dof_distribution->myp()+1), dof_distribution->np());
	}
	sprintf(fn, "%s%s", fname, pn);
  
  ifstream f(fn);
  if ((bool)f == false)
  {
    loaded_solution = false;
    return false;
  }
  
  int size;
  double v;
  
  f >> size;
  if (size != q.size())
  {
    log(W_PROG, "Solution file %s is incompatible with problem size.\n", fname);
    loaded_solution = false;
    return false;
  }
  
  for (int i=0; i<size; f >> v, q(i) = v, i++);
  
  f.close();
  
  loaded_solution = true;
  return true;
}

void CStateProblem::saveSolution(const char *fname)
{
  FILE *f;
  char fn[255], pn[255];


  
  if (fname == 0)
  {
    if (params.find("savebase") == params.end())
    {
      if (tau>0)
      {
        log(W_SOLVE, "Solution will not be saved.\n");
      }
      else
      {
        log(W_PROG, "Solution will not be saved.\n");
      }
      return;
    }
    else
    {
    	  if (dof_distribution->np() == 1)
    	  {
    		  sprintf(pn, "%s", "");
    	  }
    	  else
    	  {
    		  sprintf(pn,  ".p%dof%d", (dof_distribution->myp()+1), dof_distribution->np());
    	  }

    	  sprintf(fn, "%s.sol%s", params["savebase"].c_str(), pn);
    	  fname = fn;
    }
  }
  
  f = fopen(fname, "w");
  
  fprintf(f, "%d\n", q.size());
  int offset = 0;
  for (int c=0; c<(int)U.size(); c++)
  {
    for (int i=offset; i<offset+U[c]->ip3d->getMeshNdof(); i++)
    {
      fprintf(f, "%.15e ", q(i).getVal());
    }
    fprintf(f, "\n");
    offset += U[c]->ip3d->getMeshNdof();
  }
  
  fclose(f);
  
  if (tau>0)
  {
    log(W_SOLVE, "Solution saved to %s.\n", fname);
  }
  else
  {
    log(W_PROG, "Solution saved to %s.\n", fname);
  }
}

void CStateProblem::saveGnuplotSolution(const char *fname)
{
  char fn[255];
  
  calculateSolution();
  
  if (fname == 0)
  {
    if (params.find("savebase") == params.end())
    {
      if (tau>0)
      {
        log(W_SOLVE, "Solution will not be saved.\n");
      }
      else
      {
        log(W_PROG, "Solution will not be saved.\n");
      }
      return;
    }
    else
    {
      sprintf(fn, "%s.gpt", params["savebase"].c_str());
      fname = fn;
    }
  }
  
  output->saveGnuplot(fname);
  
  if (tau>0)
  {
    log(W_SOLVE, "Solution saved to %s.\n", fname);
  }
  else
  {
    log(W_PROG, "Solution saved to %s.\n", fname);
  }
}

void CStateProblem::saveEdgeSolution(const char *edge_name)
{
  const char *fname = (params["savebase"] + "_" + edge_name + ".dat").c_str();
  ofstream f(fname);
  int eid = mesh->getPhysNameId(edge_name);
  
  int N = output->size();
  
  calculateSolution();

  f << "# Description of output:\n"
    << "# ";
  for (int j=0; j<dim; j++) f << "x" << j+1 << " ";
  for (unsigned int c=0; c<Uinfo.size(); c++) f << Uinfo[c]->name.c_str() << " ";
  f << endl;
  
  for (int e=0; e<mesh->getNedges(); e++)
  {
    if (mesh->getEdgePhysicalNo(e) == eid)
    {
      for (int i=0; i<mesh->getNedgeDofs(); i++)
      {
	int ind = mesh->getEdgeNode(e, i);
	for (int j=0; j<dim; j++) f << mesh->getCoord(ind, j).getVal() << " ";
	ind = mesh->getNvolumes()*mesh->getNvolumeDofs() + mesh->getNfaces()*mesh->getNfaceDofs() + e*mesh->getNedgeDofs() + i;
	for (unsigned int c=0; c<U.size(); c++)
	{
	  f << U[c]->node_value[ind] << " ";
	  for (int j=0; j<dim; j++) f << U[c]->node_grad[N*j+ind] << " ";
	}
	f << endl;
      }
    }
  }
  
  f.close();
  log(W_PROG, "Solution at edge '%s' saved to %s\n", edge_name, fname);
}

void CStateProblem::log(int window, const char *format, ...)
{
  static char buffer[1024];
  va_list l;
  int rank;
  
  va_start(l, format);
  vsprintf(buffer, format, l);
  va_end(l);
  
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

  if (wm != 0 && rank == 0) wm->log(window, buffer);
  
  // print the same text to log file
  if (logfile != 0 && window == W_PROG)
  {
    fprintf(logfile, "%s", buffer);
    fflush(logfile);
  }
}

void CStateProblem::logPct(int window, int pct)
{
  if (wm != 0) wm->logPct(window, pct);
  
  if (pct >= 100)
  {
    // print the same text to log file
    if (logfile != 0 && window == W_PROG)
    {
      fprintf(logfile, "\n");
    }
  }
}


